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The  human  malaria  parasite  Plasmodium  falciparum  goes  through  a  complex  life  cycle,  including  a  roughly 
48-hour-long  intraerythrocytic  developmental  cycle  (I DC)  in  human  red  blood  cells.  A  better  understanding 
of  the  metabolic  processes  required  during  the  asexual  blood-stage  reproduction  will  enhance  our  basic 
knowledge  of  P.  falciparum  and  help  identify  critical  metabolic  reactions  and  pathways  associated  with 
blood-stage  malaria.  We  developed  a  metabolic  network  model  that  mechanistically  links  time-dependent 
gene  expression,  metabolism,  and  stage-specific  growth,  allowing  us  to  predict  the  metabolic  fluxes,  the 
biomass  production  rates,  and  the  timing  of  production  of  the  different  biomass  components  during  the  I  DC. 
We  predicted  time-  and  stage-specific  production  of  precursors  and  macromolecules  for  P.  falciparum 
(strain  HB3),  allowing  us  to  link  specific  metabolites  to  specific  physiological  functions.  For  example,  we 
hypothesized  that  coenzyme  A  might  be  involved  in  late- 1  DC  DNA  replication  and  cell  division.  Moreover, 
the  predicted  ATP  metabolism  indicated  that  energy  was  mainly  produced  from  glycolysis  and  utilized  for 
non-metabolic  processes.  Finally,  we  used  the  model  to  classify  the  entire  tricarboxylic  acid  cycle  into 
segments,  each  with  a  distinct  function,  such  as  superoxide  detoxification,  glutamate/glutamine  processing, 
and  metabolism  of  fumarate  as  a  byproduct  of  purine  biosynthesis.  By  capturing  the  normal  metabolic 
and  growth  progression  in  P.  falciparum  during  the  I  DC,  our  model  provides  a  starting  point  for  further 
elucidation  of  strain-specific  metabolic  activity,  host-parasite  interactions,  stress-induced  metabolic 
responses,  and  metabolic  responses  to  antimalarial  drugs  and  drug  candidates. 


Introduction 

Malaria  constitutes  a  major  human  health  threat,  with  more  than 
207  million  clinical  cases  and  627  000  deaths  annually.^  Despite 
the  success  of  current  anti-malarial  efforts,  including  vector  control, 
the  improvement  of  diagnostic  testing,  and  drug  development,  the 
complete  control  and  elimination  of  this  disease  is  hindered  by 
an  emerging  resistance  to  existing  drugs  and  lack  of  effective 
prophylactic  vaccines/  Thus,  efforts  to  develop  truly  novel  treat¬ 
ment  options  for  this  disease  warrant  basic  research  into  the 
different  mechanisms  used  by  Plasmodium  falciparum,  the  most 
virulent  causative  agent  of  malaria,  to  survive  and  proliferate  in 
host-based  physiological  environments/ 
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One  important  feature  of  this  process  is  the  organism’s  life 
cycle,  during  which  it  adopts  distinct  and  radically  different 
morphological  stages/  The  malaria  parasite  invades  the  host 
through  the  bite  of  an  infected  mosquito.  The  infective  sporozoite 
rapidly  moves  to  the  liver,  where  it  initially  proliferates  into 
merozoites,  largely  asymptomatically  to  the  host.  Subsequently, 
the  merozoites  enter  the  bloodstream.  While  some  merozoites 
develop  sexual  forms  that  re-infect  mosquitoes,  others  begin  a 
roughly  48  hour-long  intraeiythrocytic  developmental  cycle  (IDC) 
of  asexual  reproduction.  This  allows  the  parasites  to  infect  many 
more  red  blood  cells  and  leads  to  the  well-known  malarial 
symptom  of  recurring  fever.  During  this  cycle,  the  parasites  are 
susceptible  to  drug  treatments,  and  efforts  to  elucidate  essential 
biological  activities  such  as  metabolism,  the  focus  of  this  study, 
provide  a  foundation  for  developing  novel  therapies. 

Experimental  studies  of  P.  falciparum  metabolism  can  be 
classified  into  two  categories  as  focused  either  on  single  pathways 
or  on  more  comprehensive  “omics”-based  studies.  In  the  former 
category,  studies  have  shed  light  on  how  P.  falciparum  during 
the  IDC  satisfies  its  energy  demand  by  oxidizing  the  bulk  of  the 
taken-up  glucose  into  lactate  through  the  glycolysis  pathway.^ 
Although  the  organism  has  all  the  enzymes  to  execute  a  fully 
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functional  tricarboxylic  acid  (TCA)  cycle,  its  role  is  not  fully 
elucidated,  as  only  a  fraction  of  the  glucose  is  processed  in  this 
pathway.^’^  However,  these  studies  could  not  comprehensively 
capture  all  aspects  of  P.  falciparum  metabolism  during  the 
IDC.  Conversely,  the  latter  studies  include  high-throughput 
“omics”  profiling,  c.g.,  genome-scale  transcriptomics  profiling 
in  P.  falciparum  during  different  stages  of  its  life  cycle/’^  Of 
particular  interest  for  IDC  studies,  the  expression  data  for  nearly 
all  genes  (including  metabolic  genes)  have  been  collected  at  each 
hour  during  the  IDC,  providing  a  continual  transcriptomics 
readout  of  stage-specific  developments.^’^^  Again,  these  data 
have  primarily  been  analyzed  in  terms  of  individual  metabolic 
pathways  and  not  from  a  fully  integrated  metabolic  analysis  that 
takes  into  account  the  entire  IDC/^’^^  Here,  we  performed  a 
systems-level  investigation  of  IDC  metabolism  using  a  combi¬ 
nation  of  the  high-throughput  gene  expression  data  with  in  silica 
metabolic  network  modeling  to  systematically  connect  altered 
genetic  transcriptions  to  enzymes  and  metabolic  activities. 


systems  biology  representations  of  metabolism  use  the  consti¬ 
tuent  components  (genes,  metabolites,  and  reactions)  to  predict 
metabolic  activity  and  growth  phenotypes  of  organisms.^^  Genome- 
scale  metabolic  networks  consist  of  interconnected  biochemical 
reactions,  each  processing  particular  metabolites  spontaneously  or 
through  enzyme(s)  encoded  by  gene(s),  which,  when  they  were 
analyzed  under  certain  constrained  conditions,  such  as  limited 
nutrient  uptake,  can  predict  cellular  growth  (biomass  accumulation) 
and  other  phenotypic  functions  related  to  metabolism.  For 
instance,  metabolic  networks  for  P.  falciparum  have  been  developed 
and  used  to  identify  essential  genes/reactions  that  represent  candi¬ 
dates  for  target-based  anti-malarial  drug  discovery.  Importantly, 
these  networks  were  integrated  with  gene  expression  data  at 
different  stages  of  the  IDC,  generating  static  models  of  stage- 
specific  P.  falciparum  metabolic  activities  and  correctly  predicting 
metabolite  exchanges  between  the  parasite  and  the  host.^^’^^ 
A  limitation  with  this  approach  is  that  these  models  did  not  capture 
the  parasite’s  stage-specific  growth,  failing  to  link  genotypic 
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Fig.  1  Schematic  description  of  integrating  a  metabolic  network  with  time-series  gene  expression  data.  We  constructed  a  set  of  metabolic  states  over  a 
time  course  based  on  time-series  gene  expression  data  by  altering  nominal  fluxes  obtained  from  a  metabolic  network  that  represent  average  and  typical 
fluxes  through  each  reaction.  The  example  network  contains  five  metabolites  (A-E),  two  uptake  reactions,  three  enzymatic  reactions,  and  one  biomass 
reaction.  In  Step  I,  we  obtained  the  set  of  nominal  fluxes  that  satisfy  the  mass  balance  of  each  metabolite  and  the  typical  biomass  composition.  In  Step  II, 
we  mapped  the  time-series  gene  expression  data  to  their  corresponding  reactions.  In  Step  III,  given  the  nominal  fluxes  and  gene  expression  data,  we 
calculated  a  set  of  time-series  fluxes.  For  each  reaction,  the  time-dependent  pattern  of  these  fluxes  followed  that  of  the  corresponding  gene  expression 
data,  while  the  average  of  the  fluxes  was  equal  to  (or  as  close  as  possible  to)  the  nominal  flux.  In  Step  IV,  we  determined  the  time-series  biomass 
production  rate,  by  adding  the  production  rates  of  all  biomass  components.  We  finally  constructed  time-series  metabolic  states  that  showed  time- 
dependent  alterations  of  reaction  fluxes,  biomass  production,  and  net  production  of  each  biomass  component. 
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alterations  and  phenotypic  growth  of  P.  falciparum  during  the 
IDC.  These  models  used  fixed  biomass  functions,  and,  thus, 
were  not  able  to  capture  the  distinct  time-dependent  synthesis 
of  each  biomass  component  and  relate  these  components  to 
stage-specific  cellular  activities. 

To  overcome  these  limitations,  we  created  an  integrated  “global” 
P.  falciparum  metabolic  network  model  that  can  predict  metabolic 
fluxes,  biomass  production  rates  (representative  of  specific  growth 
rates),^^  and  the  production  of  biomass  components  at  each  hour 
during  the  IDC.  We  validated  our  model  by  comparing  the  predicted 
production  of  the  biomass  and  its  components  with  the  experimen¬ 
tally  observed  stage-specific  growth  and  macromolecular  syntheses. 
In  addition,  we  analyzed  the  predicted  net  production  of  small 
metabolites,  the  ATP  production  and  consumption,  and  the  fluxes 
through  the  TCA  cycle,  providing  a  rationale  for  the  stage-specific 
functions  of  these  molecules,  energy  metabolism,  and  metabolic 
pathways.  Thus,  our  model  provides  a  framework  to  link  gene 
expression,  metabolic  fluxes,  and  growth  phenotypes  that  can  be 
used  to  model  and  interpret  P.  falciparum  metabolism  under 
different  conditions,  including  drug  treatments  and  other 
physiological  stress  conditions. 

Results  and  discussion 

Metabolic  model  of  P.  falciparum  during  the  IDC 

We  created  a  model  that  describes  all  metabolic  activities  in 
P.  falciparum  HB3  during  asexual  reproduction  in  red  blood  cells, 
by  estimating  (l)  the  distribution  of  metabolite  flows  (or  fluxes),  (2) 
an  overall  biomass  production  rate,  and  (3)  the  production  rate  of 
each  component  of  the  organism’s  biomass  at  each  hour  of  the  IDC. 
We  used  a  modified  metabolic  network  reconstruction^^  (see  Experi¬ 
mental  section  for  details)  combined  with  hourly  gene  expression 
data  from  the  highest  time-resolution  experiment  available  to  date.^ 
This  approach  used  the  time-series  gene  expression  data^  to  alter 
a  set  of  “nominal”  metabolic  fluxes,  representing  the  average 
physiological  reaction  fluxes  of  P.  falciparum  in  the  IDC. 

Despite  the  lack  of  a  reliable  correlation  between  gene  transcrip¬ 
tion  levels  and  enzyme  activities,  a  number  of  methodologies  have 
demonstrated  the  ability  to  capture  condition-specific  metabolic 
behaviors  via  the  integration  of  a  metabolic  network  with  either 
absolute  gene  transcriptional  data^^“^^  or  differential  gene  expres¬ 
sions  compared  to  particular  reference  conditions.^’^^  Here,  we 
constructed  a  hybrid  version  of  these  approaches,  as  the  gene 
expression  data  represent  gene  transcriptional  differences  between 
the  sample  at  a  particular  time  point  and  a  mixed  pool  of  samples 
from  all  time  points.^  Importantiy,  we  used  a  semi-continuous 
approach  to  account  for  gradual  changes  in  metabolite  fluxes  and 
growth  phenotypes  rather  than  considering  metabolic  reactions  in  a 
binaiy  manner  (active  or  inactive  depending  on  predefined  stage- 
specific  metabolic  states).^^’^^  This  approach  represents  a  data- 
driven  methodology  in  which  the  time-dependent  state  of  the 
metabolic  network  (and  organism)  is  driven  and  determined  by 
alterations  in  gene  transcription. 

Fig.  1  schematically  outlines  the  major  steps  we  used  to  perform 
the  integration  of  the  metabolic  network  with  time-series  gene 


expression  data.  In  this  exemplar,  we  used  a  set  of  time-series  gene 
expression  data  and  a  metabolic  network  containing  five  meta¬ 
bolites  (A-E),  two  uptake  reactions,  three  enzymatic  reactions,  and 
one  biomass  reaction.  In  Step  I,  we  initially  used  the  network  to 
construct  a  nominal  flux  distribution  that  satisfied  the  mass 
balance  of  each  metabolite  and  the  average  biomass  composition. 
In  Step  II,  we  mapped  the  gene  expression  data  into  reactions. 
Given  the  nominal  fluxes  and  relative  expression  data,  we  generated 
a  set  of  time-series  metabolic  fluxes  that  fluctuated  around  their 
nominal  value  (Step  III).  For  each  reaction,  the  time-dependent 
pattern  of  these  fluxes  followed,  as  closely  as  possible,  that  of  the 
corresponding  expression  levels.  Using  these  fluxes,  we  could 
determine  the  syntheses  and  net  production  of  all  biomass  com¬ 
ponents  (B,  C,  and  E  in  Fig.  1),  which  when  summed  up  yielded  the 
overall  time-series  biomass  production  rate  (Step  IV). 

The  malaria  model  also  incorporated  mass  conservation 
constraints  to  account  for  the  growth  from  one  organism  to 
16-32  organisms  through  four  to  five  cell  division  cycles  during  the 
IDC.  Furthermore,  given  that  there  is  no  one-to-one  mapping 
between  gene  expression  levels  and  enzyme  activity,  we  introduced 
time  shifts  in  the  model  to  account  for  the  time  difference 
between  when  a  gene  is  transcribed  (the  available  data)  to  when 
the  synthesized  proteins  appear  and  are  active  in  the  organism 
(metabolic  activity  as  discussed  in  the  paper).  This  procedure 
resulted  in  time-series  metabolic  states  that  indicated  time- 
dependent  alterations  of  reaction  fluxes,  the  biomass  produc¬ 
tion,  and  the  synthesis  of  each  biomass  component.  We  provide 
a  detailed  description  of  these  steps  and  their  implementation  in 
the  Experimental  section. 

Prediction  of  the  overall  biomass  production 

Using  the  developed  model,  we  predicted  the  overall  biomass 
production  rate  /i  of  P.  falciparum  at  each  of  the  48  hours  during 


Fig.  2  The  predicted  overall  biomass  production  rates  fi  of  Plasmodium 
falciparum  during  the  intraerythrocytic  developmental  cycle.  The  whole 
intraerythrocytic  developmental  cycle  was  classified  into  ring,  trophozoite,  and 
schizont  stages.^°  The  unit  of  the  biomass  production  rate  is  g  h“^  gDW“^ 
denoting  gram  biomass  per  hour  per  gram  dry  weight  of  the  P.  falciparum 
merozoite  at  the  beginning  of  the  cycle. 
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the  IDC  (see  Fig.  2).  The  prediction  shows  an  overall  qualitative 
agreement  with  the  experimentally  observed  growth  pheno¬ 
types  of  the  parasite  during  the  three  stages  (ring,  trophozoite, 
and  schizont)  throughout  the  cycle.  We  predicted  that  the 
biomass  production  rate  started  from  a  veiy  low  level  and 
increased  during  the  ring  stage  (first  18  hours).  The  relatively 
low  production  rate  is  consistent  with  the  experimental  slow- 
growth  phenotype  in  the  ring  stage,^^  during  which  the  primary 
activity  of  the  organism  is  remodeling  its  internal  structures.^^ 
Further,  we  predicted  sustained  biomass  production  rates  at 
relatively  high  levels  during  the  trophozoite  stage  (18-30  hours), 
during  which  P.  falciparum  is  observed  to  grow  rapidly  by  con¬ 
suming  most  of  the  cytoplasm  of  its  host  red  blood  cells.^^  Finally, 
the  biomass  production  rates  after  30  hours  decreased,  which  is 
compatible  with  the  observed  schizont-stage  P.  falciparum  shift  in 
focus  from  growth  to  cell  division.^^  In  particular,  the  predicted 
biomass  production  rate  at  the  end  of  the  IDC  was  close  to  that  at 
the  beginning,  suggesting  that  the  metabolic  program  goes  back  to 
its  original  state  after  one  round  of  the  cycle  and  readies  itself  for 
another  round  of  infection.^^ 


Prediction  of  time-dependent  macromolecule  biomass- 
component  synthesis 

P.  falciparum  accumulates  biomass  by  taking  up  nutrients  and 
synthesizing  required  components,  such  as  DNA,  RNA,  protein, 
and  phospholipids.  However,  these  processes  are  not  uniformly 
active  during  the  IDC  and  instead  exhibit  time-  or  stage-specific 
dependencies  indicative  of  tightly  controlled  metabolic  processes 
during  merozoite  asexual  reproduction.  We  used  our  model  to 
predict  the  net  amounts  of  each  of  these  macromolecules  produced 
during  different  time  intervals  throughout  the  IDC,  and  compared 
the  predicted  results  with  the  corresponding  HB3  strain-specific 
experimental  data.^^’^^  This  comparison  was  especially  relevant  as 
the  measured  levels  from  the  experimental  system  of  P.  falciparum- 
infected  eiythrocytes  are  not  obscured  by  difficulties  in  separating 
host  and  parasite  components.  Mature  human  red  blood  cells 
contain  no  nucleus  to  synthesize  DNA  and  RNA,^^  are  thought  not 
to  synthesize  protein,^^  and  have  negligible  phospholipid  syntheses 
and  metabolism.^^  Hence,  the  measured  syntheses  of  these  mole¬ 
cules  are  attributable  to  P.  falciparum  and  not  to  the  eiythrocytes. 


A:  RNA  B:  Protein 


Time  (hour) 


D:  Phospholipids 


Simulation  Experiment 


Fig.  3  Macromolecule  syntheses  in  Plasmodium  falciparum  during  the  intraerythrocytic  developmental  cycle.  The  synthesized  amount  of  RNA  (A), 
protein  (B),  DNA  (C),  and  phospholipids  (D)  during  a  given  time  interval  was  equal  to  the  amount  of  material  synthesized  during  the  interval  normalized  by 
the  maximum  amount  produced  during  the  entire  cycle.  The  predicted  amounts  were  compared  with  the  corresponding  experimental  data  for  DNA  and 
RNA,^^  and  for  phospholipids.^^  The  horizontal  bars  indicate  the  length  of  the  time  intervals.  The  colors  of  these  bars  represent  the  simulation  results 
(blue)  and  experimental  data  (green). 
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Fig.  3A  and  B  show  that  the  predicted  syntheses  of  RNA  and 
protein,  respectively,  peaked  during  the  mid-IDC.  Specifically,  both 
the  prediction  and  experimental  data  showed  that  the  net  amount 
of  RNA  synthesized  was  relatively  high  between  12  and  36  hours 
after  infection,  with  the  highest  amount  produced  at  around 
24  hours.^^  Fig.  3B  shows  that  our  prediction  results  for  the 
synthesized  amount  of  protein  were  relatively  high  between 
18  and  40  hours  and  at  the  highest  level  between  24  and  36  hours. 
Although  no  quantitative  experimental  measurements  of  protein 
synthesis  were  available,  our  predictions  were  qualitatively 
supported  by  the  experimental  data  from  polyacrylamide  gel 
electrophoresis,  which  showed  that  the  synthesized  amounts  of 
soluble  proteins  were  relatively  high  between  24  and  36  hours, 
while  those  for  antigens  were  high  between  18  and  40  hours.^^ 
Fig.  3C  and  D  show  that  our  predictions  and  experimental 
results  indicate  that  DNA  and  phospholipids,  respectively,  were 
mainly  synthesized  in  P.  falciparum  during  the  late  IDC.  Fig.  3C 
shows  that  the  predicted  DNA  synthesis  amount  was  the  highest 
between  36  and  40  hours,  compared  to  the  maximum  experimental 


amounts,  which  occurred  a  few  hours  earlier  during  the  30-  to 
36-hour  time  interval.  Fig.  3D  shows  that  both  the  predicted 
and  experimental  synthesized  amounts  of  phospholipids 
during  the  30-  to  43-hour  time  intervals  were  higher  compared 
to  the  earlier  intervals.  Although  the  predicted  and  experimental 
amounts^^  reached  their  maximums  at  different  time  intervals,  the 
trends  in  both  sets  of  data  showed  general  increases  from  the  early 
to  the  late  IDC,  confirming  that  phospholipids  were  mainly  synthe¬ 
sized  during  the  late  stage  of  the  cycle.  The  timing  of  the  DNA  and 
phospholipid  syntheses  prepares  the  organism  for  cell  division 
during  the  schizont  stage  by  generating  the  needed  genomic 
materials  and  cellular  membranes  to  support  cell  replication. 

Despite  the  qualitative  consistency  in  the  timing  of  macro- 
molecular  production,  the  predicted  and  experimental  results  still 
exhibited  quantitative  discrepancies.  For  example,  the  predicted 
maximum  in  DNA  synthesis  occurred  six  hours  later  than  experi¬ 
mentally  observed.  These  discrepancies  are  partly  a  reflection  of  our 
limited  consideration  of  post-transcriptional  and  post-translational 
regulations.^^  As  outlined  in  the  Experimental  section,  we  did 
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Fig.  4  The  predicted  time-dependent  production  of  biomass  metabolites  of  Plasmodium  falciparum.  The  heat  map  denotes  the  predicted  time- 
dependent  production  levels  of  each  biomass  metabolite  of  P.  falciparum,  in  which  orange,  grey,  and  blue  colors  represent  high,  normal,  and  low 
production  levels,  respectively.  Based  on  the  time-dependent  production,  we  classified  these  metabolites  into  four  groups.  Groups  I,  II,  and  III  include  the 
metabolites  mainly  produced  during  the  early  (ring  stage),  middle  (trophozoite  and  early  schizont  stages),  and  late  (schizont  stage)  periods  of  the 
intraerythrocytic  developmental  cycle,  respectively,  while  Group  IV  includes  the  metabolites  for  which  the  production  levels  were  basically  constant 
throughout  the  intraerythrocytic  developmental  cycle.  The  vertical  dashed  lines  indicate  the  boundaries  between  ring,  trophozoite,  and  schizont  stages. 
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approximately  account  for  these  processes  by  shifting  the 
gene  expression  data  by  the  average  time  delay  derived  from 
comparing  a  set  of  bi-hourly  time-series  proteomics  data  with  the 
corresponding  transcriptional  levels.^^  However,  time  delays  might 
vaiy  non-linearly  during  the  developmental  stages,  e.g.,  due  to 
synchronization  switches  or  just  a  lack  of  available  ribosomes.  For 
ribosomal  proteins,  their  transcription  levels  peaked  during  the 
ring  and  early  trophozoite  stages,^  whereas  the  corresponding 
proteins  were  abundant  only  during  the  trophozoite  and  early 
schizont  stages.^^  Thus,  the  variation  in  the  number  of  ribosomal 
proteins  may  also  contribute  to  changes  in  protein  synthesis  rates 
and  in  the  time  delay  of  the  proteomics  data  compared  to  the 
corresponding  transcriptomics  data  throughout  the  IDC. 

Prediction  of  time-dependent  precursor  synthesis 

In  addition  to  the  production  of  macromolecules,  P.  falciparum 
synthesizes  a  number  of  critical  small-molecule  precursors 
which  function  as  cofactors  in  different  biological  processes. 
Fig.  4  shows  the  results  of  our  time-dependent  prediction  of  the 
precursor  synthesis  classified  into  four  general  groups.  Groups  I, 
II,  and  III  include  the  metabolites  mainly  produced  during  the 
early  (ring  stage),  middle  (trophozoite  and  early  schizont  stages), 
and  late  (schizont  stage)  periods  of  the  IDC,  respectively,  whereas 
Group  IV  includes  the  metabolites  for  which  the  production/ 
uptake  rate  was  roughly  constant  throughout  the  IDC. 

Group  I.  The  early  IDC  synthesis  included  pyridoxal  5-phosphate, 
the  active  form  of  vitamin  B6,  and  NAD,  a  ubiquitous  redox 
intermediate.  The  timing  of  these  syntheses  suggests  that  they  play 
a  key  biological  role  during  this  period,  although  both  molecules  are 
known  to  be  important  in  a  variety  of  biological  processes.  Pyridoxal 
5-phosphate  is  a  cofactor  for  more  than  100  known  enzymes,^^ 
including  13  P.  falciparum  enzymes  suggested  as  potential  druggable 
targets  by  Kronenberger  et  alf^  Given  the  predicted  high  production 
levels  of  pyridoxal  5-phosphate  in  the  ring  stage,  we  examined  the 
stage-specific  expressions  of  11  out  of  the  13  enzymes  for  which  the 
related  gene  expression  data  were  available.^  We  could  categorize  all 
enzymes  into  stage-specific  activities,  except  for  lysine  decarbo^^^lase 
(PFD0285c),  which  peaked  during  26-36  hours,  ie.,  overlapping  with 
the  late  trophozoite  and  early  schizont  stages.  Thus,  we  found  one, 
five,  and  four  enzymes  linked  to  the  ring,  trophozoite,  and  schizont 
stage,  respectively.  Table  SI  (ESIt)  provides  information  for  each 
enzyme.  In  particular,  we  could  stage-match  the  production  of 
pyridoxal  5-phosphate  with  a  cysteine  desulfurase  (encoded  by 
MAL7P1.150  [iscS])  at  the  late  ring  stage.  This  suggests  an  initial 
functional  link  between  pyridoxal  5-phosphate  and  the  iron- 
sulfur  complex  in  the  mitochondria,  which  potentially  affects 
electron  transfer,  catalysis,  and  regulatory  processes.^^ 

Group  II.  Concomitant  with  the  RNA  and  protein  syntheses 
during  the  trophozoite  and  early  schizont  stages,  a  number  of 
precursors  were  also  preferentially  produced  during  this  time, 
suggesting  a  role  for  them  in  protein  production  or  other  activities 
during  the  mid-IDG.  For  example,  de  ^ovo-synthesized  heme  was 
demonstrated  to  be  a  requirement  for  optimal  protein  syntheses.^^’^^ 
Fig.  4  also  shows  that  the  polyamines  (spermidine  and  putrescine) 
were  predicted  to  achieve  peak  production  levels  during  mid-IDC, 
yet  they  are  involved  in  a  wide  range  of  functions,  including  the 


stabilization  of  DNA,  RNA,  and  proteins.^^  However,  our  prediction 
of  their  preferential  synthesis  during  mid-IDC  is  supported  by  the 
experimental  observation  that  inhibition  of  polyamine  syntheses 
had  no  visible  effect  during  the  ring  stage  but  caused  morphological 
growth  arrest  in  the  trophozoite  stage.^^  Finally,  our  prediction 
of  the  time-dependent  synthesis  of  1 0-formyl tetrahydrofolate 
coincided  with  protein  synthesis,  tied  to  its  function  of  donating 
its  formyl  moiety  to  produce  formylmethionine,  the  initiator  of 
mRNA  translation.^^ 

Group  in.  Fig.  4  shows  that  the  peak  production  of  the  precursors 
included  in  this  group  roughly  coincided  with  the  enhanced 
production  of  DNA  and  phospholipids  during  the  schizont 
stage.  In  this  group,  the  predicted  synthesis  of  ubiquinone  was 
confirmed  by  noting  that  the  highest  experimental  ubiquinone 
concentration  occurred  in  the  schizont  stage.^^ 

The  observed  time  dependency  among  this  group  of  mole¬ 
cules  also  allowed  us  to  hypothesize  about  their  biological 
role.  For  example,  *S-adenosyl  L-methionine  involves  both  the 
trophozoite-stage  polyamine  syntheses^^  and  DNA  methylation, 
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Fig.  5  Predicted  energy  production  and  consumption.  (A)  Schematic 
description  of  energy  production  and  consumption.  Energy  (in  form  of 
ATP)  was  produced  from  glycolysis  and  other  metabolic  pathways,  and 
consumed  by  non-glycolytic  metabolism  and  non-metabolic  activity. 
(B)  Predicted  time-dependent  ATP  production  and  consumption  with  respect 
to  metabolic  and  non-metabolic  processes  (excluding  ATP  used  for  RNA 
synthesis).  The  unit  of  the  production  or  consumption  is  mmol  h“^  gDW“\ 
denoting  millimole  per  hour  per  gram  dry  weight  of  the  P.  falciparum 
merozoite  at  the  beginning  of  the  intraerythrocytic  developmental  cycle. 
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an  important  step  in  cell  division  during  the  schizont  stage.^^ 
Based  on  our  predictions  of  an  increased  production  of  5-adenosyl 
L-methionine  during  the  late  IDC,  we  hypothesized  that  the  primaiy 
function  of  this  molecule  might  be  involved  in  DNA  methylation. 
A  further  example  was  coenzyme  A,  whose  inclusion  in  Group  III 
suggests  that  its  function  is  related  to  schizont-stage  activities. 
Coenzyme  A  is  reported  to  be  necessaiy  for  DNA  replication  and 
cell  division  in  fission  yeast, but  its  hypothesized  role  in  DNA 
replication  of  P.  falciparum  remains  to  be  confirmed.  Conversely, 
coenzyme  A  biosynthesis  is  a  valid  antimicrobial  drug  target  and, 
based  on  its  essentiality  for  P.  falciparum  in  red  blood  cells 
and  different  biochemical  characteristics  from  its  counterpart 
in  humans,^^  may  be  considered  as  a  potential  anti-malarial 
drug-development  target. 

Group  IV.  This  group  includes  molecules  for  which  we 
predicted  roughly  constant  production/uptake  rates  throughout 
the  IDC.  Because  the  production/uptake  rates  were  driven  by 
altered  gene  transcription  profiles,  our  model  did  not  capture 
other  time-  and  stage-specific  processes.  Similarly,  any  lack  of  a 
direct  gene-reaction  relationship  in  our  model  would  prevent 


us  from  modeling  such  a  reaction.  For  example,  in  the  network 
model,  riboflavin  was  directly  taken  up  by  the  parasite  and  used 
for  biomass  production.^^  However,  as  no  gene  was  associated 
with  this  process,  we  did  not  model  any  time-dependent  varia¬ 
tion  in  riboflavin  accumulation. 

Prediction  of  energy  production  and  consumption 

It  has  long  been  proposed  that  P.  falciparum  obtains  energy  by 
anaerobically  metabolizing  glucose  into  lactate  through 
the  glycolysis  pathway.^  However,  it  remains  unclear  whether 
other  metabolic  pathways  contribute  appreciably  to  energy 
production  and  what  biological  processes  use  the  produced 
energy.  Using  our  model,  we  calculated  the  time-dependent 
ATP  production  from  glycolysis,  the  production  and  consump¬ 
tion  by  other  metabolic  pathways,  and  the  consumption  by 
non-metabolic  activities  (z.e.,  conversion  of  ATP  to  ADP).  Fig.  5 
shows  that  compared  to  glycolysis,  other  metabolic  pathways 
produced  negligible  amounts  of  ATP,  and  that  these  non-glycolytic 
pathways  consumed  much  less  ATP  than  non-metabolic  activities. 
In  addition,  these  ATP  production  and  consumption  pathways 


Fig.  6  Predicted  metabolic  fluxes  through  the  tricarboxylic  acid  cycle.  (A)  Reactions  in  the  Plasmodium  falciparum  tricarboxylic  acid  (TCA)  cycle.  The 
arrows  represent  the  directions  of  the  predicted  fluxes  and  different  line  styles  indicated  our  classification  of  the  TCA  cycle  into  several  segments,  each  of 
which  had  a  distinct  function.  (B)  The  predicted  time-dependent  flux  profiles  of  these  reactions  during  the  intraerythrocytic  developmental  cycle  (IDC). 
(C)  The  rescaled  flux  profile  of  aconitate  hydratase  (ACONT),  citrate  synthase  (CS),  and  isocitrate  dehydrogenase  (ICDH)  during  the  IDC.  The  unit  of  the 
fluxes  is  mmol  h“^  gDW“\  denoting  millimole  per  hour  per  gram  dry  weight  of  the  P.  falciparum  merozoite  at  the  beginning  of  the  IDC.  AKGDH,  oc- 
ketoglutarate  dehydrogenase;  FUM,  fumarase;  MDH,  i-malate  dehydrogenase;  SUCD,  succinate  dehydrogenase;  and  SUCOAS,  succinate-coenzyme  A 
ligase. 
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reached  their  highest  levels  during  the  trophozoite  and  late 
schizont  stages,  indicating  that  P.  falciparum  relied  almost 
exclusively  on  the  glycolysis  pathway  to  generate  energy 
during  the  IDC,  most  of  which  was  used  for  non-metabolic 
activities,  such  as  RNA  production,  protein  syntheses,  and 
cellular  growth. 

Prediction  and  analysis  of  metabolite  fluxes  in  the  TCA  cycle 

P.  falciparum  possesses  a  mitochondrial  compartment  and 
encodes  all  of  the  enzymes  in  the  TCA  cycle,^  indicating  the 
presence  of  a  fully  functional  TCA  cycle.  However,  the  function  of 
the  cycle  still  remains  unclear,  because  P.  falciparum  does  not 
generate  the  bulk  of  its  energy  by  processing  glucose  through  the 
TCA  cycle.^’^  Instead,  P.  falciparum  metabolizes  glucose  into 
lactate,  which  is  subsequently  secreted.^’^  Here,  we  predicted 
the  time-series  fluxes  for  all  the  reactions  in  P.  falciparum  during 
the  IDC  to  clarify  the  metabolic  role  of  the  TCA  cycle. 

The  predicted  fluxes  through  the  TCA  reactions  shared  some 
common  features.  For  example.  Fig.  6A  shows  that  all  reaction 
fluxes  were  in  the  oxidative  direction,  which  was  consistent 
with  isotope  profiling  experiments.^’^  In  addition.  Fig.  6B  and  C 
show  that  the  predicted  flux  profiles  for  all  the  reactions  in  the 
TCA  cycle  peaked  during  the  schizont  stage  (30-48  hours), 
indicating  the  importance  of  these  reactions  for  this  stage  of 
the  IDC. 

The  TCA  reaction  fluxes  also  showed  different  time-dependent 
profiles,  suggesting  that  we  could  identify  individual  segments 
with  different  biological  functions.  Fig.  6B  shows  that  the 
fluxes  through  the  reactions  catalyzed  by  citrate  synthase 
(CS),  aconitate  hydratase  (ACONT),  and  isocitrate  dehydro¬ 
genase  (ICDH)  were  much  lower  than  other  reactions.  This 
segment  of  the  TCA  cycle  involves  the  detoxification  of  super¬ 
oxide  leaked  from  the  mitochondrial  respiration  chain  where 
the  leaked  superoxide  is  converted  by  dismutase  into  hydrogen 
peroxide,  which  is  subsequently  reduced  to  water.^^  This 
reduction  of  hydrogen  peroxide  oxidizes  NADPH  into  NADP, 
which,  in  turn,  is  reduced  by  ICDH.^^’^^  Furthermore,  the 
similarity  in  the  flux  profiles  of  a-ketoglutarate  dehydrogenase, 
succinate-coenzyme  A  ligase,  and  succinate  dehydrogenase 
(AKGDH,  SUCOAS,  and  SUCD  in  Fig.  6,  respectively)  enabled 
us  to  group  these  reactions  into  one  segment,  which  we  linked 
to  oxidation  of  2-oxoglutarate  generated  from  glutamate  and 
glutamine.^  This  oxidation  is  continued  by  fumarase  (FUM  in 
Fig.  6),^  which  also  metabolizes  fumarate  from  the  purine 
pathway.^^  Finally,  we  predicted  that  part  of  the  produced 
L-malate  was  transported  out  of  mitochondria  while  the 
remaining  L-malate  was  converted  in  the  TCA  cycle  into  oxalo- 
acetate  by  L-malate  dehydrogenase  (MDH  in  Fig.  6). 

Conclusion 

The  IDC  represents  the  major  growth-through-replication  phase 
of  malaria  parasites  in  humans,  and  its  roughly  48-hour  cycle 
corresponds  to  the  well-known  recurring  host-fever  attacks. 
During  this  phase,  the  merozoite  undergoes  three  distinct 


transformations  that  allow  it  to  absorb  nutrients  from  the 
infected  blood  cell  and  undergoes  four  or  five  cycles  of  replica¬ 
tion  to  produce  16  or  32  new  merozoites,  respectively.  To 
examine  how  P.  falciparum  executes  and  adjusts  its  metabolism 
during  this  replication  phase,  we  created  a  genome-scale  meta¬ 
bolic  network  model  that  predicted  metabolic  fluxes,  the  bio¬ 
mass  production  rate,  and  the  production  of  each  biomass 
metabolite  at  each  hour  during  the  IDC.  This  work  builds  on 
and  extends  the  capabilities  of  the  existing  metabolic  models  of 
blood-stage  P.  falciparum}^'^^  In  particular,  we  have  extended 
the  metabolic  network  models  to  (l)  predict  time-dependent 
biomass  production  and  (2)  allow  for  a  time-variable  biomass 
function,  so  that  we  could  correctly  model  and  capture 
the  stage-specific  synthesis  of  each  biomass  metabolite  and 
capture  stage-specific  growth  phenotypes.  Our  results  indicate 
that  the  model  provides  detailed  mechanistic  links  between 
transcriptomics,  metabolism,  and  phenotypically  distinct 
organism  development.  This  lays  the  foundation  for  the  future 
interpretation  and  exploration  of  P.  falciparum  metabolism 
under  a  wide  variety  of  stress  conditions. 

Experimental 

Metabolic  network  of  P.  falciparum 

From  the  two  recently  developed  metabolic  networks  of 
P.  falciparum,^^’^^  we  started  our  model  development  using 
the  /TH366  gene-protein-reaction  formulations,  stoichiometric 
representations  of  all  reactions,  and  biomass  objective  functions.^^ 
Although  equivalent  in  scope,  the  zTH366  network  differentiates 
the  reactions  catalyzed  by  isozymes  from  those  by  multiple-unit 
protein  complexes, while  the  PlasmoNet  network  assumes 
that  all  of  these  reactions  are  catalyzed  by  isozymes. The  more 
extensive  gene-protein-reaction  formulation  theoretically 
allows  for  a  finer  distinction  when  integrating  the  gene  expres¬ 
sion  data  with  the  individual  metabolic  reactions. 

Modifications.  We  made  modifications  to  the  existing  zTH366 
set  of  reactions  by  changing  the  compartment  of  glutathione 
reductase  from  mitochondria  to  cytosol,^^’^^  adding  the  genes 
encoding  phosphoethanolamine  methyltransferase,^^  and  the 
enzymes  associated  with  the  hemoglobin  degradation.^^’^^ 
We  also  added  a  new  set  of  reactions,  including  serine 
decarbo^^^lase,^^  methylenetetrahydrofolate  reductase,^^  cardio- 
lipin  synthase,^^  the  polymerization  of  hemoglobin-generated 
heme,  the  degradation  of  the  heme  by  glutathione,  mito¬ 
chondrial  thioredoxin  reductase,  and  mitochondrial  peroxire- 
doxin.^^’^^’^^  We  only  allowed  uptake  and  secretion  for  the  set  of 
metabolites  that  were  found  to  be  transported  by  P.  falciparum 
during  the  IDC  (see  the  complete  list  in  Table  S2  in 
We  also  assumed  free  secretion  of  CO2  and  free  transports  of  H2O 
and  H^,  and  we  added  sulfate  uptake  and  urea  secretion  to  avoid  a 
zero  biomass  production  rate. 

We  redefined  the  original  biomass  objective  function  into  a 
new  set  of  biomass  functions  (denoted  by  B)  to  allow  for 
monitoring  biomass  metabolite  production  at  different  times 
during  the  IDC.  Thus,  the  new  formulation  B  of  the  biomass 
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function  describes  the  biological  activities  necessary  for  the 
biomass  production  in  P.  falciparum,  through  the  following 
equations: 

AA«  +  tRNAa  AA-tRNA«,  a  =  1,. .  .,20  (l) 

^  CAAitRNAi  +  •  •  •  +  CAA20tRNA20  +  Protein  (2) 

Cy^TpATP  “I"  Cf^'ppCTP  "I"  Cq'ppGTP  "I"  Cu^ppUXP 

— ^  RNA  "I"  (c^^pp  "I"  C(]"pp  "I"  Cq"pp  "I"  Cu^p^Diphosphate 

(3) 

CdATP^^'T'P  ^  ^dCTpdC'T'P  "t  CdGTpdGr'T'P  "t"  C^TTP^ -I'-l'P 

DNA  +  (CdATP  ^dcTP  ^dGTP  CdTTp)Diphosphate 

(4) 

Cp^PC  "I"  CppPE  "t"  CpgPS  "I"  CpjPI  "I"  Cspg;SPH  "I"  CglCL  — ^  f  (^) 

^  (j)  (for  each  metabolite  i  that  was  a  polyamine, 
a  cofactor,  or  an  inorganic  ion)  (6) 

^energyH20  +  CgnergyATP 

^  ^^energyADP  +  CgnergyPhoSphatC  +  CenergyH^  (?) 

Eqn  (1)  and  (2)  represent  the  addition  of  amino  acids  for 
protein  syntheses,  where  AA  denotes  an  amino  acid,  a  repre¬ 
sents  its  index,  tRNA^  denotes  tRNA  used  to  transport  amino 
acid  a,  AA-tRNA^  indicates  the  complex  of  amino  acid  a  and  its 
corresponding  tRNA,  and  Caa^  is  the  corresponding  stoichio¬ 
metric  coefficient  for  amino  acid  a.  Similarly,  eqn  (3)  and  (4) 
represent  the  use  of  nucleotides  for  the  syntheses  of  RNA 
and  DNA,  respectively.  Eqn  (5)  accounts  for  the  addition  of 
the  phospholipids  phosphatidylcholine,  phosphatidylethanol- 
amine,  phosphatidylserine,  phosphatidylinositol,  sphingomye¬ 
lin,  and  cardiolipin  (denoted  as  PC,  PE,  PS,  PI,  SPH,  and  CL  in 
eqn  (5),  respectively),  which  are  all  necessary  components  of 
cellular  membranes.^^’^^  Eqn  (6)  and  (7)  denote  the  utilization 
of  other  metabolites  (M^)  and  energy  consumption  for  the 
biomass  production,  respectively.  In  the  above  equations, 
c  represents  the  coefficient  of  a  metabolite  in  the  original 
biomass  objective  function.  One  should  note  that  carbo¬ 
hydrates  are  not  part  of  the  original  biomass  objective  function, 
primarily  because  there  are  no  experimental  data  available. 
This  omission  should  not  be  a  major  limitation  for  the  meta¬ 
bolic  network  formulation  for  blood-stage  malaria,  as  the 
parasites  do  not  store  energy  reserves  during  this  stage,  but 
instead  rely  heavily  on  glucose  uptake  (C.  Plata,  personal 
communication).  Once  experimental  data  are  available  for 
other  metabolites,  such  as  carbohydrates,  one  could  add  their 
compositions  into  the  biomass  objective  function. 

We  have  provided  the  constructed  model  representing 
P.  falciparum  metabolism  during  the  IDC  (in  MATLAB  format) 
in  the  Protocol  SI  in  ESI.t 

Gene  expression  data  for  P.  falciparum  during  the  IDC 

To  capture  time-dependent  metabolism  as  accurately  as 
possible,  we  used  the  time-series  gene  expression  data  at  the 


highest  time  resolution,  collected  hourly  from  synchronized 
populations  of  P.  falciparum  HB3  during  the  IDC.^’^o  data 
set  is  measured  using  two-channel  microarrays  in  which  the 
treatment  channel  contained  the  sample  at  one  time  point, 
while  the  control  channel  contained  a  mixed  pool  of  samples 
from  all  time  points. Basing  our  model  construction 
and  calculations  on  data  for  the  HB3  strain^  enabled  us  to 
validate  our  results  using  experimentally  measured  synthesis 
rates  of  DNA,  RNA,  and  protein  that  are  only  available  for 
this  strain.^^ 

Calculation  of  metabolic  fluxes  in  P.  falciparum  at  each  hour 
during  the  IDC 

Fig.  1  shows  the  overall  scheme  for  calculating  metabolic  fluxes 
of  P.  falciparum  at  each  time  point  of  the  IDC.  We  used  a 
stepwise  approach  to  construct  a  set  of  time-series  reaction 
fluxes  that  fluctuated  around  their  nominal  values  based  on 
corresponding  time-series  gene  expression  data. 

Step  I.  We  calculated  a  set  of  nominal  fluxes  that  were  equal 
to  the  most  parsimonious  metabolite  flow  through  the  meta¬ 
bolic  network,  based  on  the  assumption  of  the  most  efficient 
utilization  of  nutrients.^^  This  was  done  by  solving  two  mini¬ 
mization  problems,  the  first  of  which  was  the  minimization  of 
the  overall  nutrient  uptake  rate, 

^vy/a; 

j€T 


lb  <  V  <  ub 

Vj  =  for  each  j  e  B, 

where  v  denotes  the  flux  vector  with  element  Vj  indicating 
the  flux  through  reaction  j  in  units  of  mmol  h“^  gDW“^, 
or  millimole  per  hour  per  gram  dry  weight  of  the  original 
merozoite  (the  asexual  form  of  P.  falciparum  capable  of  infect¬ 
ing  erythrocytes);  T  represents  the  set  of  nutrient  uptake  reac¬ 
tions;  ay  denotes  the  human  serum  concentration^^  of  the 
nutrient  taken  up  by  reaction  j;  S  indicates  the  stoichiometric 
matrix;  lb  and  ub  denote  the  lower  and  upper  bounds  of  the 
fluxes,  respectively;  and  represents  the  average  growth  rate 
of  P.  falciparum  during  the  IDC. 

The  coefficients  (l/ay)  in  the  objective  function  ensure  a 
preference  to  uptake  nutrients  associated  with  high  concentra¬ 
tions.  We  set  the  average  growth  rate  based  on  the  observation 
that  during  the  IDC,  the  organism  undergoes  an  average  of  four 
to  five  cell  divisions.^®  Because  the  total  number  of  divisions  is 
not  controlled  under  experimental  conditions,  the  total  mass  at 
the  end  of  the  IDC  is  |  x  (2^  +  2^)  or  24  times  the  initial 
merozoite  mass.  This  sets  the  value  for  p^  to  (24  —  l)/48  or 
0.48  gram  new  biomass  per  hour  per  gram  dry  weight  of  the 
original  merozoite  weight  (g  h“^  gDW“^).  By  setting  the  con¬ 
straints  Vj  =  p^  (for  j  e  B),  we  fixed  a  certain  flux  through  each 
biomass  function,  which  ensures  that  the  mass  of  each  bio¬ 
mass  component  at  the  end  of  the  IDC  is  24  times  its  initial 
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mass,  corresponding  to  the  observation  that,  on  average,  one 
merozoite  generates  24  copies  of  itself  during  the  IDC. 

Given  minimum  nutrient  uptake,  we  subsequently  minimized 
the  sum  of  intracellular  fluxes  by  solving  the  following  problem: 

min  ^|v,- 

MT 

s.t.  S  •  V  =  0 
lb  <  V  <  ub 

Vj  =  for  each  reaction  j  e  B 
for  each  reaction  i  G  T, 

where  v/  represents  the  minimum  nutrient  uptake  rates 
obtained  in  eqn  (8),  and  the  inequalities  Vj  <  Vj*  (for  j  e  T) 
indicate  that  the  nutrient  uptakes  are  constrained  by  the 
corresponding  minimum  values  v/. 

After  the  two  minimizations,  we  obtained  a  set  of  nominal 
fluxes  representative  of  typical  physiological  fluxes  through 
the  reactions.  For  example,  fluxes  through  the  glycolysis  pathway 
(including  lactate  production  and  secretion)  were  among  the 
largest,  consistent  with  the  well-established  fermentative  glucose 
utilization  by  the  parasite.^ 

Step  II.  Next,  we  processed  the  time-series  gene  expression 
data  of  P.  falciparum  during  the  IDC^  to  obtain  an  expression  level 
Vj  for  each  reaction  j  at  each  time  point  t.  After  first  smoothing  the 
original  gene  expression  data  using  the  locally  weighted  scatterplot 
smoothing  (LOWESS)  method  with  a  span  parameter  of  30%,^ 
we  subsequently  mapped  the  resultant  gene  expression  data 
to  reactions  by  established  methods.^^’^^  If  a  single  enzyme 
catalyzed  a  reaction,  we  assigned  the  corresponding  expression 
level  of  that  gene  to  that  reaction.  If  a  reaction  was  catalyzed  by 
multiple  isozymes,  we  assigned  their  maximum  gene  expres¬ 
sion  level  to  that  reaction.  If  a  reaction  was  catalyzed  by  a 
protein  complex  composed  of  multiple  subunits,  we  assigned 
their  minimum  gene  expression  level  to  that  reaction. 

Time  shift  between  gene  transcription  and  protein  translation. 
In  this  article,  we  have  discussed  metabolism  in  terms  of  when 
enzymatic  reactions  occur  and  when  the  model  predicts  synthesis 
of  biomaterials.  As  the  available  experimental  data  of  gene 
transcription  levels  do  not  directly  give  us  protein  levels,  we 
time-shifted  the  available  transcription  levels  by  a  time  difference 
equal  to  the  difference  between  peak  expression  levels  and  peak 
protein  abundance  levels.^^  For  reactions  associated  with  directly 
measured  time-dependent  protein  abundances,  the  time  delays 
were  taken  as  the  differences  between  the  peak  times  of  protein 
abundances  and  those  of  the  corresponding  mRNA  levels.  If  these 
reactions  were  part  of  the  same  metabolic  pathway  and  proteomic 
data  were  unavailable  for  some  enzymes  in  the  same  pathway,  the 
time  delay  for  the  reactions  with  missing  data  was  assigned  based 
on  the  average  peak-time  differences  of  the  reactions  with  experi¬ 
mentally  available  data.  For  all  remaining  reactions,  the  time 
delays  were  11  hours,  corresponding  to  the  median  of  peak  time 
differences  between  all  the  measured  protein  abundances  and  the 
corresponding  transcriptomics  data.^^ 


Linear  normalization  of  expression  data.  Finally,  we  normalized 
the  expression  data  linearly,  so  that  the  mean  and  minimum  of 
each  reaction’s  time-series  expression  values  were  one  and  zero, 
respectively.  This  assumes  that  the  minimum  expression  value 
corresponded  to  complete  inactivation  of  the  reaction,  and  the 
mean  value  corresponded  to  the  nominal  flux  through  the  reaction. 

Step  III.  Given  the  nominal  fluxes  (from  Step  I)  and  the 
expression  level  rj  for  each  reaction  j  at  each  time  point  t  (from 
Step  II),  we  calculated  the  reaction  fluxes  at  each  time  point  by  first 
minimizing  the  difference  between  the  reaction  fluxes  and  the 
product  of  their  nominal  fluxes  and  expression  values  as  follows: 

min  ^  I  vj  —  rj  •  I 

j^G 


lb  <  <  ub, 

where  vj  represents  the  flux  through  reaction  j  at  time  point  t,  and 
G  represents  the  set  of  intracellular  irreversible  reactions  that  can 
be  associated  with  gene  expression  data.  Because  the  expression 
data  of  transport  (uptake  or  secretion)  reactions  are  mapped  to 
transport  fluxes  for  molecules  involved  in  both  metabolic  (subject 
to  constraints  in  our  model)  and  non-metabolic  processes  (not 
subject  to  constraints),  G  does  not  include  these  reactions. 
Examples  of  transport  reactions  for  non-metabolic  processes 
include  water  for  regulation  of  osmotic  pressure^^’^^  and 
transport  for  maintaining  intracellular  pH  levels.^^’^^  G  also 
excludes  reactions  that  can  carry  fluxes  in  both  directions, 
because  we  cannot  unambiguously  determine  the  directions  of 
their  fluxes.  In  total,  we  excluded  eight  such  reactions  among  the 
317  intracellular  reactions  that  were  associated  with  gene  expres¬ 
sion  data  and  not  in  dead-end  pathways. 

Given  that  the  problem  formulated  in  eqn  (10)  may  have 
multiple  solutions,  we  selected  the  solution  closest  to  the  nominal 
flux  distribution,  based  on  the  assumption  that  reaction  fluxes 
minimally  deviated  from  their  nominal  values.  We  did  this  by 
solving  the  following  additional  optimization  problem: 

min  ^  I  vj  —  vj^  | 

J 


s.t.  S  V  =  0 

lb  <  <  ub 


(11) 


J^G 

where  Jf  is  the  optimal  value  for  the  objective  function  from 
the  previous  optimization  problem  defined  by  eqn  (10),  and  the 
last  constraint  ensures  that  this  solution  is  one  of  the  optimal 
solutions  for  eqn  (10). 

After  solving  eqn  (10)  and  (11),  we  obtained  metabolic  flux 
distributions  for  all  time  points  during  the  IDC.  These  time-series 
fluxes  varied  around  their  nominal  values  in  time-dependent 
patterns  similar  to  those  of  the  corresponding  expression  data. 
If  no  gene  expression  data  were  available,  we  assigned  the  reaction 
the  time-independent  value  of  the  nominal  flux  distribution. 
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The  intent  of  the  above  approach  was  to  be  able  to  use  the 
relative  gene  expression  level  as  a  semi-quantitative  indicator 
that  the  cell  is  readying  itself  to  transcribe  the  produced  mRNA  into 
proteins  to  execute  or  affect  some  biological/metabolic  process. 
Thus,  the  metabolic  model  captures  time-dependent  aspects  of 
when  the  gene  transcripts  became  available,  as  well  as  the  relative 
magnitudes  of  the  expression  levels.  We  further  accounted  for  the 
mRNA/protein  mismatch  by  using  the  experimentally  measured 
time  difference  between  when  a  gene  was  transcribed  and  when 
the  synthesized  proteins  appeared  in  the  organism.  The  use  of 
relative  expression  levels  around  the  average  levels  allows  the 
optimization  procedure  flexibility  in  determining  metabolic  fluxes 
consistent  with,  but  not  equal  to,  gene  expression  levels. 

Step  IV.  We  determined  the  overall  biomass  production  level 
fi  at  each  time  point  t  from 


M  =■ 


E  vjwy 

Ew,- 

J^B 


(12) 


where  Wj  indicates  the  biomass  fraction  of  the  metabolite(s) 
associated  with  the  biomass  function  j.  We  defined  Wj  as 


Wj  = 


CijWi 

^  1000’ 

I 


(13) 


where  Cy  represents  the  coefficient  of  metabolite  i  in  biomass 
function  j,  Wt  denotes  the  molecular  weight  of  the  metabolite, 
and  the  factor  1000  converts  mol  into  mmol. 


Simulation  environment 

We  constructed  the  model  and  ran  the  simulations  in  MATLAB 
(2012a,  MathWorks,  Natick,  MA)  using  the  COBRA  toolbox.^^  The 
metabolic  model  of  P.  falciparum  during  the  IDC  (in  MATLAB 
format)  is  provided  as  Protocol  SI  in  ESI.t 
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